其他
读取GEO数据库的单细胞转录组表达矩阵文本文件的一种方式
最近在读AUCell包的文档,链接是:http://bioconductor.org/packages/release/bioc/html/AUCell.html,这个包的教程我已经写完了, 在 :使用AUCell包的AUCell_calcAUC函数计算每个细胞的每个基因集的活性程度
发现AUCell包使用了 GSE60361 数据集的单细胞转录组表达矩阵,是直接读取文本文件文件,代码具有学习价值,值得反复分享,如下:
library(GEOquery)
attemptsLeft <- 20
while(attemptsLeft>0)
{
geoFile <- tryCatch(getGEOSuppFiles("GSE60361", makeDirectory=FALSE), error=identity)
if(methods::is(geoFile,"error"))
{
attemptsLeft <- attemptsLeft-1
Sys.sleep(5)
}
else
attemptsLeft <- 0
}
gzFile <- grep(".txt.gz", basename(rownames(geoFile)), fixed=TRUE, value=TRUE)
message(gzFile)
txtFile <- gsub(".gz", "", gzFile, fixed=TRUE)
message(txtFile)
gunzip(filename=gzFile, destname=txtFile, remove=TRUE)
library(data.table)
geoData <- fread(txtFile, sep="\t")
geneNames <- unname(unlist(geoData[,1, with=FALSE]))
exprMatrix <- as.matrix(geoData[,-1, with=FALSE])
rm(geoData)
dim(exprMatrix)
rownames(exprMatrix) <- geneNames
exprMatrix[1:5,1:4]
# Remove file
file.remove(txtFile)
# Save for future use
mouseBrainExprMatrix <- exprMatrix
save(mouseBrainExprMatrix, file="exprMatrix_AUCellVignette_MouseBrain.RData")
当然啦, 每个人写代码风格不一样,我就不会这样写。
最后读入的表达矩阵被整理好了,是小鼠的约2万个基因的3千多个细胞的表达矩阵,如下所示:
> dim(exprMatrix)
[1] 19972 3005
>
> exprMatrix[1:5,1:4]
1772071015_C02 1772071017_G12 1772071017_A05 1772071014_B06
Tspan12 0 0 0 3
Tshz1 3 1 0 2
Fnbp1l 3 1 6 4
Adamts15 0 0 0 0
Cldn12 1 1 1 0
总有人问GEO表达矩阵如何下载如何读取
自己多看一些优秀代码,这个很具有学习价值的代码分享给大家哈!其实初学者遇到的绝大部分问题仅仅是因为基础知识不扎实而已。再怎么强调生物信息学数据分析学习过程的计算机基础知识的打磨都不为过,我把它粗略的分成基于R语言的统计可视化,以及基于Linux的NGS数据处理:
把R的知识点路线图搞定,如下:
了解常量和变量概念 加减乘除等运算(计算器) 多种数据类型(数值,字符,逻辑,因子) 多种数据结构(向量,矩阵,数组,数据框,列表) 文件读取和写出 简单统计可视化 无限量函数学习
Linux的6个阶段也跨越过去 ,一般来说,每个阶段都需要至少一天以上的学习:
第1阶段:把linux系统玩得跟Windows或者MacOS那样的桌面操作系统一样顺畅,主要目的就是去可视化,熟悉黑白命令行界面,可以仅仅以键盘交互模式完成常规文件夹及文件管理工作。 第2阶段:做到文本文件的表格化处理,类似于以键盘交互模式完成Excel表格的排序、计数、筛选、去冗余,查找,切割,替换,合并,补齐,熟练掌握awk,sed,grep这文本处理的三驾马车。 第3阶段:元字符,通配符及shell中的各种扩展,从此linux操作不再神秘! 第4阶段:高级目录管理:软硬链接,绝对路径和相对路径,环境变量。 第5阶段:任务提交及批处理,脚本编写解放你的双手。 第6阶段:软件安装及conda管理,让linux系统实用性放飞自我。